R语言中dnorm, pnorm, qnorm与rnorm以及随机数

发表于 | 分类于
| 字数统计: 3,642 | 阅读时长 ≈ 17

前言

在R语言中，与正态分布（或者说其它分布）有关的函数有四个，分别为dnorm,pnorm,qnorm和rnorm，其中，dnorm表示密度函数，pnorm表示分布函数，qnorm表示分位数函数，rnorm表示生成随机数的函数。在R中与之类似的函数还有很多，具体的可以通过help(Distributions)命令去查看，对于分位数或百分位数的一些介绍可以看这篇笔记《分位数及其应用》，关于正态分布的知识可以看这篇笔记《正态分布笔记》。

现在这篇笔记就介绍一下这些函数的区别。

R中的随机数背景

R提供了多种随机数生成器(random number generators, RNG)，默认采用的是Mersenne twister方法产生的随机数，该方法是由Makoto Matsumoto和Takuji Nishimura于1997年提出来的，其循环周期是$2^{19937}-1$。R里面还提供了了Wichmann-Hill、Marsaglia-Multicarry、Super-Duper、Knuth-TAOCP-2002、Knuth-TAOCP和L’Ecuyer-CMRG等几种随机数生成方法，可以通过RNGkind()函数进行更改，例如，如果要改为WIchmann-Hill方法，就使用如下语句：

1
RNGkind(kind="Wich")

set.seed()随机数种子

在R中使用随机数函数，例如rnorm()函数来生成的随机数是不一样的，有时我们在做模拟时，为了比较不同的方法，就需要生成的随机数都一样，即重复生成相同的随机数，此时就可以使用set.seed()来设置随机数种子，其参数为整数，如下所示：

1
2
3
4
5
6
7
8
9
> set.seed(1)
> runif(5)
[1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819
> set.seed(2)
> runif(5)
[1] 0.1848823 0.7023740 0.5733263 0.1680519 0.9438393
> set.seed(1)
> runif(5)
[1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819

dnorm

dnorm中的d表示density，norm表示正态贫，这个函数是正态分布的概率密度(probability density)函数。

正态分布的公式如下所示：

给定x，μ和σ后，dnorm()这个函数返回的就是会返回上面的这个公式的值，这个值就是Z-score，如果是标准正态分布，那么上述的公式就变成了这个样子，如下所示：

现在看一个案例，如下所示：

1
2
3
dnorm(0,mean=0,sd=1)
# 这个是标准正态分布函数
# [1] 0.3989423

dnorm(0,mean=0,sd=1)由于是标准正态分布函数的概率密度，这个命令其实可以直接写为dnorm(0)即可，如下所示：

1
2
dnorm(0)
# [1] 0.3989423

再看一个非标准正态分布的案例，如下所示：

1
2
dnorm(2, mean = 5, sd = 3)
# [1] 0.08065691

虽然在dnorm()中，x是一个概率密度函数(PDF,Probability Density Function)的独立变量，但它也能看作是一组经过Z转换后的一组变量，现在我们看一下使用dnorm来绘制一个正态分布的概率密度函数曲线，如下所示：

1
2
3
4
5
6
7
8
z_scores <- seq(-3,3,0.1)
# 生成一个Z-score的向量
z_scores
# [1] -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5
# [17] -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1  0.0  0.1
# [33]  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.0  1.1  1.2  1.3  1.4  1.5  1.6  1.7
# [49]  1.8  1.9  2.0  2.1  2.2  2.3  2.4  2.5  2.6  2.7  2.8  2.9  3.0

现在使用dnorm()函数计算一下Z_scores的概率密度，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
dvalues <- dnorm(z_scores)
dvalues
# [1] 0.004431848 0.005952532 0.007915452 0.010420935 0.013582969 0.017528300
# [7] 0.022394530 0.028327038 0.035474593 0.043983596 0.053990967 0.065615815
# [13] 0.078950158 0.094049077 0.110920835 0.129517596 0.149727466 0.171368592
# [19] 0.194186055 0.217852177 0.241970725 0.266085250 0.289691553 0.312253933
# [25] 0.333224603 0.352065327 0.368270140 0.381387815 0.391042694 0.396952547
# [31] 0.398942280 0.396952547 0.391042694 0.381387815 0.368270140 0.352065327
# [37] 0.333224603 0.312253933 0.289691553 0.266085250 0.241970725 0.217852177
# [43] 0.194186055 0.171368592 0.149727466 0.129517596 0.110920835 0.094049077
# [49] 0.078950158 0.065615815 0.053990967 0.043983596 0.035474593 0.028327038
# [55] 0.022394530 0.017528300 0.013582969 0.010420935 0.007915452 0.005952532
# [61] 0.004431848

现在绘图，如下所示：

1
2
3
4
plot(z_scores,dvalues, # Plot where y = values and x = index of the value in the vector
     type = "l", # Make it a line plot
     main = "pdf of the Standard Normal",
     xlab= "Z-score")

从上面的结果可以看出，在每个Z-score处，dnorm可以绘制出这个Z-score对应的正态分布的pdf的高度。

pnorm

pnorm函数中的p表示Probability，它的功能是，在正态分布的PDF曲线上，返回从负无穷到q的积分，其中这个q指的是一个Z-score。现在我们大概就可以猜测出pnorm(0)的值是0.5，因为在标准正态分布曲线上，当Z-score等于0时，这个点正好在标准正态分布曲线的正中间，那么从负无穷到0之间的曲线面积就是整个标准正态分布曲线下面积的一半，如下所示：

1
2
3
pnorm(0)
# pnorm()默认的参数与dnorm()一样，都是标准正态分布，即平均数为0，标准差为1的正态分布
# [1] 0.5

pnorm函数还能使用lower.tail参数，如果lower.tail设置为FALSE，那么pnorm()函数返回的积分就是从q到正无穷区间的PDF下的曲线面积，因此我们就知道了，pnorm(q)与1-pnorm(q,lower.tail=FALSE)的结果是一样的，如下所示：

1
2
3
4
5
6
7
8
9
10
11
pnorm(2)
# [1] 0.9772499
pnorm(2, mean = 5, sd = 3)
# [1] 0.1586553
pnorm(2, mean = 5, sd = 3, lower.tail = FALSE)
# [1] 0.8413447
1 - pnorm(2, mean = 5, sd = 3, lower.tail = FALSE)
# [1] 0.1586553

在计算机出现之前的时代里，统计学家们使用正态分布进行统计时，通常是要查正态分布表的，但是，在计算机时代，通常都不使用正态分布表了，在R中，pnorm()这个函数完全可以取代正态分布表了，现在我们使用一个Z-scores的向量来计算一下相应的累积概率，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
# 此处还使用前面生成的z_scores
z_scores
# [1] -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5
# [17] -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1  0.0  0.1
# [33]  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.0  1.1  1.2  1.3  1.4  1.5  1.6  1.7
# [49]  1.8  1.9  2.0  2.1  2.2  2.3  2.4  2.5  2.6  2.7  2.8  2.9  3.0
pvalues <- pnorm(z_scores)
pvalues
# [1] 0.001349898 0.001865813 0.002555130 0.003466974 0.004661188 0.006209665
# [7] 0.008197536 0.010724110 0.013903448 0.017864421 0.022750132 0.028716560
# [13] 0.035930319 0.044565463 0.054799292 0.066807201 0.080756659 0.096800485
# [19] 0.115069670 0.135666061 0.158655254 0.184060125 0.211855399 0.241963652
# [25] 0.274253118 0.308537539 0.344578258 0.382088578 0.420740291 0.460172163
# [31] 0.500000000 0.539827837 0.579259709 0.617911422 0.655421742 0.691462461
# [37] 0.725746882 0.758036348 0.788144601 0.815939875 0.841344746 0.864333939
# [43] 0.884930330 0.903199515 0.919243341 0.933192799 0.945200708 0.955434537
# [49] 0.964069681 0.971283440 0.977249868 0.982135579 0.986096552 0.989275890
# [55] 0.991802464 0.993790335 0.995338812 0.996533026 0.997444870 0.998134187
# [61] 0.998650102
plot(pvalues, 
     xaxt = "n", # 不绘制x轴的刻度
     type = "l", # 使用曲线将各点连接起来
     main = "标准正态分布的CDF曲线",
     xlab= "分位数(Quantiles)",
     ylab="概率密度(Probability Density)") 
# 以下是添加x轴的刻度
# These commands label the x-axis
axis(1, at=which(pvalues == pnorm(-2)), labels=round(pnorm(-2), 2))
axis(1, at=which(pvalues == pnorm(-1)), labels=round(pnorm(-1), 2))
axis(1, at=which(pvalues == pnorm(0)), labels=c(.5))
axis(1, at=which(pvalues == pnorm(1)), labels=round(pnorm(1), 2))
axis(1, at=which(pvalues == pnorm(2)), labels=round(pnorm(2), 2))

以上就是标准正态分布的累积分布函数(CDF,Cumulative Distribution Function)曲线。

qnorm

简单来说，qnorm是正态分布累积分布函数(CDF,Cumulative Distribution Function）的反函数，也就是说它可以视为pnorm的反函数，这里的q指的是quantile，即分位数。

使用qnorm这个函数可以回答这个问题：正态分布中的第p个分位数的Z-score是多少？

现在我们来计算一下，在正态分布分布中，第50百分位数的Z-score是多少，如下所示：

1
2
3
4
5
6
7
qnorm(0.5)
# [1] 0
# 这里计算的就是在标准正态分布中，第50百分位数的Z-score是多少
pnorm(0)
# [1] 0.5
# 这里计算的就是在标准正态分布中，当Z-score是0时，它对应的曲线下面积是多少，也就是对应的是哪个百分位数

再来看一个案例：在正态分布中，第96个百分位的Z-score是多少，如下所示：

1
2
qnorm(.96)
# [1] 1.750686

再来看一个案例：在正态分布中，第99个百分位的Z-score是多少，如下所示：

1
2
qnorm(0.99)
# [1] 2.326348

再来看一下pnorm()这个函数，如下所示：

1
2
3
4
5
pnorm(qnorm(0))
# [1] 0
pnorm(2.326348)
# [1] 0.99

从上面我们可以看到，pnorm这个函数的功能是，我们知道某个Z-score是多少，它位于哪个分位数上。

接着我们进一步举例来说明一下qnorm和pnorm的具体功能，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
oldpar <- par()
par(mfrow=c(1,2))
# 设置绘图页面，页面布局是一行两列
quantiles <- seq(0, 1, by = .05)
# 以5%为步长，生成0到1的百分数
quantiles
# [1] 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75
# [17] 0.80 0.85 0.90 0.95 1.00
qvalues <- qnorm(quantiles)
# 计算每个百分位数对应的Z-score
qvalues
# [1]       -Inf -1.6448536 -1.2815516 -1.0364334 -0.8416212 -0.6744898 -0.5244005
# [8] -0.3853205 -0.2533471 -0.1256613  0.0000000  0.1256613  0.2533471  0.3853205
# [15]  0.5244005  0.6744898  0.8416212  1.0364334  1.2815516  1.6448536        Inf

现在进行绘图，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
plot(qvalues,
     type = "l", # We want a line graph
     xaxt = "n", # No x-axis
     xlab="概率密度(Probability Density)",
     ylab="Z-scores")
# Same pnorm plot from before
plot(pvalues, # Plot where y = values and x = index of the value in the vector
     xaxt = "n", # Don't label the x-axis
     type = "l", # Make it a line plot
     main = "标准正态分布的CDF曲线",
     xlab= "分位数(Quantiles)",
     ylab="概率密度(Probability Density)") 
# 绘制x轴的刻度
axis(1, at=which(pvalues == pnorm(-2)), labels=round(pnorm(-2), 2))
axis(1, at=which(pvalues == pnorm(-1)), labels=round(pnorm(-1), 2))
axis(1, at=which(pvalues == pnorm(0)), labels=c(.5))
axis(1, at=which(pvalues == pnorm(1)), labels=round(pnorm(1), 2))
axis(1, at=which(pvalues == pnorm(2)), labels=round(pnorm(2), 2))

rnorm

rnomr()函数的功能用于生成一组符合正态分布的随机数，在学习各种统计学方法时，rnorm这个函数应该是最常用的，它的参数有n,mean，sd，其中n表示生成的随机数，mean与sd分别表示正态分布的均值与标准差，现在举个例子，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
set.seed(1000)
# 设定随身数种子
rnorm(5)
# 生成5个服从标准正态分布的随机数
# [1] -0.44577826 -1.20585657  0.04112631  0.63938841 -0.78655436
n10 <- rnorm(10, mean = 70, sd = 5);n10
# 生成10个，服从均值为70，标准差为5的正态分布的随机数
# [1] 71.80351 60.47042 73.38245 74.64034 73.20814 75.98877 78.17864 68.84888 72.23209
# [10] 81.19844
n100 <- rnorm(100, mean = 70, sd = 5);n100[1:10]
# 生成100个，服从均值为70，标准差为5的正态分布的随机数
# [1] 62.09470 72.75453 65.42062 68.65545 76.84043 74.17830 75.40543 65.32326 70.71987
# [10] 70.81957
n10000 <-  rnorm(10000, mean = 70, sd = 5);n10000[1:10]
# 生成1000个，服从均值为70，标准差为5的正态分布的随机数
# [1] 66.14695 70.16345 69.54554 76.15484 74.54789 71.78985 67.85345 73.85163 67.58083
# [10] 73.98425

现在我们绘制一下上面的几个向量的直方图，看一下它们的均值是否在70附近，如下所示：

1
2
3
4
5
6
7
oldpar <- par()
par(mfrow=c(1,3))
# breaks为设置直方图的宽度
hist(n10, breaks = 5)
hist(n100, breaks = 20)
hist(n10000, breaks = 100)

在R语言中，生成不同分布的各种类型的函数都是以d,p,q,r开头的，使用原理跟上面的正态分布都一样。

runif-生成均匀分布的随机数

1
2
3
4
runif(5)        # 生成 5 个介于 0 和 1 之间的均匀分布的随机数
runif(5, 1,10)  # 生成 5 个介于 0 和 10 之间的均匀分布的随机数
rnorm(5)        # 生成 5 个正态分布的随机数，它们的中位数为 0，标准差为 1
rnorm(5, 3, 7)  # 生成 5 个正态分布的随机数，它们的中位数为 3，标准差为 7

R语言中的随机数

sample()函数是一个用于生成随机数的重要的核心函数，如果仅传递一个数值n给它，就会返回一个从1到n的自然数的排列，如果传递是n:m就是生成从n到m的随机数，如是是7,5，则会生成5个小于7的随机数，如下所示：

1
2
3
4
5
6
> sample(7)
[1] 5 1 2 6 4 3 7
> sample(7:5)
[1] 5 6 7
> sample(7, 5)
[1] 6 4 1 3 2

从上面的结果可以看出来，这些数字都是不同的，也就是说，sample函数默认情况下是不重复抽样，每个值只出现一次，如果允许有重复抽样，需要添加参数replace = TRUE，如下所示：

1
2
> sample(7, 10, replace = TRUE)
 [1] 3 7 7 2 5 7 1 4 5 5

sample函数通常会从某些向量中随机挑一些参数，如下所示：

1
2
> sample(colors(), 5)
[1] "honeydew4"    "wheat4"       "deepskyblue3" "lightgray"    "royalblue1"

也可以挑日期，如下所示：

1
2
> sample(.leap.seconds, 4)
[1] "1994-07-01 08:00:00 CST" "1977-01-01 08:00:00 CST" "1982-07-01 08:00:00 CST" "1972-07-01 08:00:00 CST"

R语言中其它分布的随机数

分布 中文名称 R中的表达式 参数
Beta 贝塔分布 beta(a,b) shape1、shape2
Binomial 二项分布 binom(n,p) size、prob
Cauchy 柯西分布 cauchy() location、scale
Chi-square 卡方分布 chisq(df) df
Exponential 指数分布 exp(lambda) rate
F F分布 f(df1,df2) df1、df2
Gamma 伽玛分布 gamma() shape、rate
Geometric 几何分布 geom() prob
Hypergeometric 超几何分布 hyper() m、n、k
Logistic 逻辑分布 logis() location、scale
Negative binomial 负二项分布 nbinom() size、prob
Normal 正态分布 norm() mean、sd
Multivariate normal 多元正态分布 mvnorm() mean、cov
Poisson 泊松分布 pois() lambda
T t分布 t() df
Uniform 均匀分布 unif() min、max
Weibull 威布尔分布 weibull() shape、scale
Wilcoxon 威尔考可森分布 wilcox() m、n

上述分布函数前面加上r，p、q、d就可以表示相应的目的：

  • r：生成相应分布的随机数；
  • d：生成相应分布的密度函数；
  • p：生成相应分布的累积概率密度函数；
  • q：生成相应分布的分位数函数。

参考资料

  1. Introduction to dnorm, pnorm, qnorm, and rnorm for new biostatisticians
  2. R数据分析——方法与案例详解.电子工业出版社出版.方匡南 朱建平 姜叶飞.2015